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Abstract 

The  linear  stability  of  circular  Couette  flow  between  concentric  infinite  cylinders  is 
considered  for  the  case  when  the  inner  cylinder  is  rotated  at  a constant  angular  velocity 
and  the  outer  cylinder  is  driven  sinusoidally  in  time  with  zero  mean  rotation.  This 
configuration  was  studied  experimentally  by  Walsh  and  Donnelly.  The  critical  Reynolds 
numbers  calculated  from  linear  stability  theory  agree  well  with  the  experimental  values, 
except  at  large  modulation  amplitudes  and  small  frequencies.  The  theoretical  values  are 
obtained  using  Floquet  theory  implemented  in  two  distinct  approaches:  1)  a truncated 
Fourier  series  representation  in  time  and  2)  a fundamental  solution  matrix  based  on  a 
Chebyshev-pseudospectral  representation  in  space.  For  large  amplitude,  low  frequency 
modulation,  the  linear  eigenfunctions  axe  temporally  complex,  consisting  of  a quiescent 
interval  followed  by  rapid  change  in  the  perturbed  flow  velocities. 
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1.  Introduction 


In  1923  G.I.  Taylor  considered  the  stability  of  the  steady  Couette  flow  between  rotating  cylin- 
ders and  discovered  a transition  to  steady  axisymmetric  toroidal  rolls  (“Taylor  vortices”), 
obtaining  agreement  between  the  measured  values  of  the  critical  rotation  rates  and  his  the- 
oretical predictions  [1].  In  the  last  thirty  years,  a substantial  amount  of  research  has  been 
conducted  on  various  aspects  of  Taylor-Couette  flow,  including  the  effect  of  time-periodic 
forcing  of  the  cylinder  rotation  rates.  Beginning  with  the  experimental  work  of  Donnelly 
et  al  [2],  [3],  where  it  was  observed  that  temporal  modulation  of  the  inner  cylinder  angular 
velocity  provides  stabilization,  the  problem  has  served  as  an  important  model  for  gaining 
understanding  into  the  effects  of  time-periodic  forcing. 

The  present  study  is  an  outgrowth  of  a research  effort  to  understand  the  interaction  of 
a solidification  interface,  and  its  associated  interfacial  instabilities,  with  some  of  the  funda- 
mental hydrodynamic  instabilities  [4].  These  problems  are  relevant  to  the  area  of  materials 
processing,  in  particular,  growth  of  a solid  from  the  melt.  For  the  Taylor-Couette  prob- 
lem, recent  studies  by  McFadden  et  al  have  demonstrated  that  a solidification  interface  can 
significantly  alter  the  centrifugal  instability  for  steady  rotation  [5],  [6],  and  the  effect  of  time- 
periodic  rotation  is  now  being  addressed.  In  developing  the  linear  theory  for  the  interaction 
of  modulated  Taylor-Couette  flow  with  a solidification  interface,  the  isothermal  problem  sub- 
ject to  time-periodic  oscillation  was  used  to  test  the  solution  procedures.  After  reviewing  the 
literature  on  this  problem,  it  was  found  that  some  discrepancies  remain  between  the  existing 
experimental  results  and  theoretical  predictions. 

Carmi  and  Tustaniwskyj  [7]  performed  a comprehensive  study  of  the  linear  theory  onset 
conditions  for  sinusoidal  torsional  oscillation  of  the  cylinders  using  Floquet  theory.  The 
study  was  comprehensive,  in  that  several  variations  of  forcing  conditions  were  considered, 
with  and  without  mean  rotation,  for  a finite-gap  width  and  for  a gap  width  that  approximates 
the  narrow-gap  limit.  Their  analysis  predicts  destabilization  predominantly  as  a result  of 
modulation.  Donnelly  and  colleagues  have  conducted  several  of  the  more  recent  experimental 
studies  on  the  effects  of  temporal  modulation.  For  the  particular  case  of  steady  rotation  of 
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the  inner  cylinder  and  sinusoidal  torsional  oscillation  of  the  outer  cylinder  about  a zero  mean, 
Walsh  and  Donnelly  [8]  observed  stabilization  due  to  the  modulation,  while  the  linear  theory 
of  Carmi  and  Tustaniwskyj  [7]  predicts  destabilization.  A recent  theoretical  analysis  by 
Barenghi  and  Jones  [9]  yields  much  better  agreement  between  linear  theory  and  experiment 
for  one  set  of  the  Walsh  and  Donnelly  experiments;  however,  this  set  of  conditions  was  not 
the  focus  of  their  study,  and  the  experimental  conditions  yielding  the  greatest  stabilization 
were  not  considered. 

The  present  study  provides  a more  complete  comparison  with  the  cases  considered  ex- 
perimentally by  Walsh  and  Donnelly  [8].  A thorough  theoretical  investigation  is  conducted 
based  on  linear  theory  employing  Floquet  analysis  in  two  independent  approaches  to  the 
numerical  solution.  A detailed  description  of  the  linear  theory  eigenmodes,  both  spatial  and 
temporal  is  included.  Although  limited  to  linear  behavior,  the  comprehensive  study  provides 
better  understanding  into  the  complex  nature  of  the  response  to  temporally-periodic  forcing 
in  the  Taylor-Couette  problem. 

2.  Theory 

Briefly,  the  governing  equations  are  the  continuity  equation  and  the  incompressible  Navier- 
Stokes  equations  for  the  velocity,  u,  and  the  pressure,  p, 

V • u = 0,  (la) 

o 1 

~ + (u  ■ V)u  + -Vp  = 1/V2u,  (lb) 

ot  p 

where  p is  the  density  and  v is  the  kinematic  viscosity.  For  cylindrical-polar  coordinates 
(r,0,2),  the  velocity  components  are  u = ( u,v,w ) in  the  radial,  azimuthal,  and  axial  direc- 
tions, respectively.  The  fluid  occupies  the  annular  region  R\  < r < R2,  where  R\  is  the  inner 
cylinder  radius  and  R2  is  the  outer  cylinder  radius.  For  the  stability  analysis,  the  cylinder 
is  assumed  to  be  infinite  in  the  axial  direction. 

We  focus  the  analysis  on  time-periodic  rotation  of  the  outer  cylinder  about  a zero  mean, 
while  the  inner  cylinder  rotates  with  a constant  angular  velocity  The  angular  velocity 
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of  the  outer  cylinder  is  assumed  to  be  Cl? (t)  = cos  ut,  where  u>  is  the  forcing  frequency. 
The  boundary  conditions  on  the  velocity  components  are  the  no-slip  and  no-penetration 
conditions  at  solid  boundaries;  therefore,  the  azimuthal  velocity  at  the  inner  cylinder  is 
i’i  = Ri^-i  and  at  the  outer  cylinder  is  V2  = #2^2(0-  The  remaining  velocity  components 
are  zero  at  both  boundaries. 

We  next  choose  dimensionless  variables.  The  length  scale  is  chosen  to  be  the  fluid  gap 
width  d = i?2  — Ri,  the  time  scale  is  chosen  to  be  d2/^,  and  the  velocity  scale  is  chosen  to  be 
v / d.  We  retain  the  same  notation  for  all  variables,  which  will  henceforth  be  dimensionless. 
The  choice  of  scaling  introduces  the  Reynolds  number  in  the  boundary  condition,  which  is 
given  by  Re  = Q.\dP /v.  The  fluid  region  then  occupies  the  range  77/ ( 1 — 77)  < r < 1/(1  — 77), 
where  77  = R1/R2. 

3.  Base  State 

The  base  state  is  given  by  u = 0,  v = v^(r,t),  and  w = 0.  The  dimensionless  azimuthal 
velocity  satisfies 


dv(°)  _ ^dV°)  ldu<0) 


,(0) 


dt 


\ dr 2 r dr 


with  the  boundary  conditions, 


v(0){t}/(1  - r j),t)  = 

1 - 77 


and 


u(0)(l/(l  - 77), t)  = 


eRe 
1 - 77 


cos  ujt. 


(2a) 


(2b) 


(2c) 


The  solution  to  the  above  linear  problem  for  the  base  flow  is  written  as  the  superposition 
of  a steady  and  a periodic  component,  v ^ = Vs(r)  + Vp(r,t).  The  steady  solution  is  found 
by  setting  the  right-hand  side  of  Eq.  (2a)  equal  to  zero,  and  applying  the  steady  part  of  the 
boundary  conditions.  The  steady  component  of  the  base  flow  is 


V, : = 


R, 


1 


1 — 772  L (1  — 77 )2  r 


772r] 


(3) 
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An  analytic  solution  for  the  periodic  component  of  the  base  flow  solution  can  be  written  in 
terms  of  Kelvin  functions.  However,  because  the  solution  to  the  stabihty  problem  must  be 
determined  numerically,  it  is  often  more  efficient  to  calculate  the  periodic  component  of  the 
base  flow  numerically  in  the  same  manner  used  to  calculate  the  perturbed  flow. 


4.  Linearized  equations 


For  the  linear  stability  analysis  of  the  time-dependent  base  state,  the  total  flow  field  variables 
are  written  as  the  superposition  of  the  base  state  component  and  a perturbation.  The 
perturbed  quantities  are  Fourier  analyzed  in  the  azimuthal  and  axial  directions,  so  that  we 
write 

u(r,z,(f),t ) 
v{r,z,(f>,t) 
w{r,z,4>,t) 


{ 


\ 


\ p(r,*,<M)  / 


( 


0 


\ 


id0*(r,  t ) 

0 

V P(0)(M)  ) 


( u(r,t ) ^ 


+ 


w(r,  t) 

\ P(r,t)  ) 


exp  + iaz), 


(4) 


where  a is  the  nondimensional  axial  wave  number  and  n is  the  azimuthal  mode  number.  The 
base  state  components  are  denoted  by  a zero  superscript  (e.g.,  td0*)  and  the  quantities  u,  v, 
w , p are  the  perturbation  amplitudes.  Note  that  only  axisymmetric  disturbances  ( n = 0) 
will  be  used  to  obtain  the  stability  results  presented  here;  however,  the  analysis  is  formulated 
for  the  general  nonaxisymmetric  case. 

The  governing  equations  for  the  perturbation  quantities  are  obtained  by  substituting  the 
above  expansion  into  the  dimensionless  form  of  Eq.  (1)  and  linearizing  in  the  perturbation 
quantities: 


u inv 

DH-+ + iaw  = 0, 

r r 


du  . td°)  A -tv?*  Du 

— — K m u — 2 v + Dp  = D2u  -} 

ot  r r r 


n2u 


u 2 inv 

a2u , 

~2  r2  ’ 


dv 


,(°) 


-r — |-  in v + t)u  -I — = D2v  + 

ot  r r 


Dv  n2v 


dw  . A . A 2 A Dw  n2w 

— + in w + tap  = D w H — 

ot  r r r* 


— a2v + 


2 * 

— aw , 


2 inu 


(5a) 

(5b) 

(5c) 

(5d) 
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for  77/(1  — 77)  < r < 1/(1  — 77).  Here  D = djdr  and  D.  = d / dr  4-  1 jr. 

The  order  of  the  system  is  reduced  by  differentiating  Eq.  (5a)  with  respect  to  r and 
subtracting  it  from  Eq.  (5b)  yielding, 

du  . A A „ n2u  2a  ,rnDv  nv 

— — h in u — 2 u + Dp  = a u — z 1 + aDw  . 

Ot  r r rz  rr2 

The  set  of  equations  (5a),  (5c),  (5d),  and  (6)  yield  a sixth-order  system  which  is  closed  by 
the  six  boundary  conditions  u = v = w = 0 imposed  at  r = 77/(1  — 77)  and  r = 1/(1  — 77). 

Alternatively,  one  can  eliminate  w and  p from  the  set  of  equations  by  straightforward 
manipulation  to  obtain  a fourth-order  equation  in  r for  the  quantity  u and  a second-order 
equation  for  v.  This  is  a commonly-used  formulation  for  hydrodynamic  stability  problems,  so 
the  equations  are  not  given  here  (see  Carmi  and  Tustaniwskyj  [7]  for  the  full  nonaxisymmetric 
equations).  The  boundary  conditions  in  this  form  become  v = u = Du  = 0 at  the  inner  and 
outer  boundaries.  As  discussed  in  the  next  section,  two  different  numerical  approaches  are 
used  to  solve  the  stability  problem.  In  the  first  approach,  the  set  of  equations  retaining  all 
the  variables  shown  above  is  used,  while  in  the  second  approach,  the  reduced  set  of  equations 
in  terms  of  u and  v alone  is  advantageous. 

5.  Numerical  Treatment 

The  set  of  equations  for  the  perturbation  quantities  are  linear  partial  differential  equations  in 
one  space  variable  and  time.  The  equations  contain  time  periodic  coefficients,  and  so  Floquet 
theory  can  be  used  to  investigate  the  stability  of  the  system.  Two  distinct  implementations 
of  Floquet  theory  are  employed.  In  the  first  approach,  the  time-periodic  part  of  the  solution 
is  represented  by  a truncated  complex  Fourier  series  in  time.  This  approach  was  employed 
by  Hall  [10]  and  by  Seminara  and  Hall  [11]  for  similar  stability  problems.  In  this  approach, 
a perturbation  quantity  /(r,  t)  is  represented  by  the  product  of  a periodic  Fourier  series  and 
an  exponential  term  with  complex  growth  rate  cr, 

f(r,t)  = e"  Y.  /4r)e'm“‘, 

|m|<M 
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The  real  part  of  cr  determines  the  stability  of  the  base  solution.  Substitution  of  the  series 
expansion  into  the  equations  and  boundary  conditions  above  yields  a large  set  of  coupled 
two-point  boundary  value  problems  in  the  spatial  variable  r for  the  Fourier  coefficients, 
fm(r).  The  set  of  equations  is  complex  and  contains  the  growth  rate  a as  a parameter. 

The  set  of  coupled  two-point  boundary  value  problems,  subject  to  the  homogeneous 
boundary  conditions,  yields  an  eigenvalue  problem  that  is  solved  using  the  computer  code 
SUPORT  in  a similar  manner  as  the  one  described  in  [12].  The  SUPORT  code  [13]  solves  the 
set  of  linear  coupled  twTo-point  boundary  value  problems  using  superposition  of  numerically 
integrated  solutions.  An  orthonormalization  procedure  is  used  to  assure  linear  independence 
of  the  intermediate  solutions.  The  user  has  a choice  of  variable-step  integration  schemes. 
For  all  the  results  presented,  a high-order  Adams-type  procedure  was  used. 

The  base  state  subject  to  periodic  forcing  is  linearly  stable  if  for  all  disturbances  crr  < 0, 
where  crr  is  the  real  part  of  cr.  For  the  condition  of  neutral  stability  (crr  = 0),  the  solution  of 
the  eigenvalue  problem  for  the  full  set  of  coupled  complex  Fourier  modes  yields  the  Reynolds 
number  and  the  imaginary  part  of  the  growth  rate  (a,)  for  fixed  values  of  the  remaining 
parameters.  Using  the  full  set  of  equations  for  the  problem  posed  here,  cr,  is  found  to  be  either 
zero  (synchronous  response)  or  lo/2  (subharmonic  response).  By  fixing  the  value  of  cr,  for 
either  synchronous  or  subharmonic  response,  symmetry  relationships  between  the  negative 
and  positive  index  Fourier  modes  can  be  used  to  halve  the  number  of  unknowns  required 
to  solve  for  the  Reynolds  number  alone.  For  axisymmetric  disturbances,  a set  of  12A/  + 12 
ordinary  differential  equations  for  the  coupled,  complex  temporal  Fourier  modes  must  be 
solved.  The  value  of  M used  depends  on  the  parameter  values,  in  particular,  the  value  of 
the  forcing  frequency.  Lower  frequency  calculations  are  more  temporally  complicated  and 
require  a greater  number  of  Fourier  modes  in  time.  For  the  results  presented,  the  value  of  M 
ranged  from  4 to  24,  depending  on  frequency.  A discussion  of  the  accuracy  and  convergence 
of  this  approach  is  contained  in  the  next  section. 

The  periodic  part  of  the  base  flow  solution  is  obtained  by  assuming 

Vp{r,t)  = Real[V(r)e'wt } 


-7- 


where  the  real  part  corresponds  to  cosine  forcing  at  the  boundaries.  It  is  through  the 
periodic  time  dependence  of  the  base  flow  solution  terms  in  the  perturbation  equations  that 
the  individual  temporal  Fourier  solution  modes  are  coupled.  Substituting  the  above  form 
into  the  equation  and  boundary  conditions  given  by  Eq.  (2)  yields  a two-point  boundary 
value  problem  for  the  complex  amplitude  function  V(r),  which  has  an  analytic  solution  in 
the  form  of  Kelvin  functions  of  the  first  and  second  kind.  Because  a variable  stepsize  is  used 
to  integrate  the  perturbation  equations  numerically,  it  is  computationally  more  efficient  to 
solve  for  V(r)  using  SUPORT,  rather  than  evaluate  the  Kelvin  functions  which  comprise 
the  analytic  solution.  Calculating  the  complex  amplitude  function  V{r)  numerically  in  this 
manner  does  not  degrade  the  accuracy  of  the  required  base  flow  solution. 

The  second  numerical  approach  employed  to  solve  the  stability  problem  consists  of  ap- 
proximating the  spatial  behavior  of  the  solutions,  using  the  equation  set  for  u and  v alone, 
via  the  pseudospectral  technique  in  the  physical  domain,  as  described  in  references  [14]  and 
[15].  The  approach  corresponds  to  expanding  the  solutions  in  terms  of  truncated  series  of 
Chebyshev  polynomials  Tn(s), 

N 

u(r,t)  = ^2  UM  Tn(s), 

n=  0 
N 

v{r,t)  = J2  vn{t)Tn{s), 

n=  0 

where  s = 2(r  — 77/(1  — 7/))  — 1.  The  pseudospectral  discretization  requires  that  the  solution 
expansions  satisfy  the  governing  equations  at  specific  collocation  points  for  the  Chebyshev 
polynomials.  When  implemented  in  the  physical  domain,  the  unknowns  are  the  solution 
values  at  the  collocation  points,  instead  of  the  expansion  coefficients  as  in  purely  spectral 
methods. 

In  the  pseudospectral  approach,  the  spatial  differential  operators  in  the  governing  partial 
differential  equations  are  replaced  by  discrete  matrix  operators.  As  a result,  the  governing  set 
of  partial  differential  equations  and  boundary  conditions  becomes  a set  of  coupled  ordinary 
differential-algebraic  equations  in  time  for  the  unknown  solution  values  at  the  collocation 
points.  The  algebraic  equations  result  from  the  boundary  conditions.  Because  the  size  of  the 
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set  of  coupled  equations  depends  on  the  number  of  original  partial  differential  equations,  the 
size  is  minimized  by  using  the  equations  formulated  in  terms  of  u and  v alone,  as  discussed 
in  the  previous  section;  the  total  number  of  unknowns  in  this  formulation  is  2 (N  + 10- 
In  this  numerical  solution  approach,  the  computational  effort  for  determining  the  periodic 
component  of  the  base  flow  by  numerical  integration  or  by  evaluation  of  the  Kelvin  functions 
is  equivalent,  since  the  amplitude  values  V{r)  are  required  only  at  fixed  collocation  points 
(i.e.,  they  need  only  be  calculated  once  and  stored). 

The  computer  code  DASSL  [16]  is  used  to  integrate  the  differential-algebraic  system  in 
time.  The  algorithm  uses  backward  differentiation  with  up  to  fifth-order  accuracy  obtained 
using  a multi-step  approach.  The  implicit  integration  procedure  is  appropriate  for  the  inte- 
gration of  the  equation  set  obtained  from  a Chebyshev  spatial  discretization.  Also,  because 
differential-algebraic  systems  often  behave  like  systems  of  stiff  differential  equations,  they 
require  that  special  error  estimation  schemes  be  used  in  order  to  obtain  stable,  accurate 
solutions  [17]. 

In  this  second  approach,  Floquet  analysis  is  implemented  by  constructing  a K x I\ 
fundamental  solution  matrix,  where  K is  the  number  of  differential  equations  (A'  = 2 N — 4). 
The  K columns  of  this  matrix  are  linearly  independent  calculated  solutions  for  the  unknowns 
at  the  end  of  one  forcing  period.  The  eigenvalues  of  this  matrix  are  the  Floquet  multipliers 
from  which  the  complex  growth  rate  a is  obtained.  The  advantage  of  the  second  approach 
is  that  it  yields  K values  of  a resulting  from  the  discrete  solution  modes,  providing  several 
eigenvalues  that  are  accurate  approximations  to  the  spectrum.  Since  periodic  forcing  can 
excite  different  temporal  response  (i.e.,  subharmonic  versus  synchronous),  knowledge  of  the 
relative  stability  of  various  modes  provided  by  the  a values  simplifies  the  determination  of 
the  critical  stability  boundaries  in  parameter  space. 

The  two  distinct  numerical  solution  approaches  provide  an  independent  check  on  the  cal- 
culated results.  In  addition,  each  of  the  two  approaches  has  advantages  and  disadvantages 
depending  on  the  nature  of  the  solution.  For  example,  if  the  spatial  variation  is  limited  more 
to  the  boundary  regions  of  the  domain,  the  SUPORT  approach  with  its  variable  mesh  incre- 


-9- 


ment  can  optimize  the  computation  effort  for  a given  global  accuracy;  while  for  complicated 
temporal  behavior,  the  time  integration  scheme  of  the  second  approach  may  be  better,  since 
the  temporal  Fourier  representation  may  converge  slowly  and  require  an  impractical  number 
of  terms  for  reasonable  accuracy. 

6.  Results 

As  discussed  in  the  introduction,  the  focus  of  the  present  study  is  a comprehensive  com- 
parison of  linear  stability  theory  predictions  for  the  onset  conditions  with  the  experimental 
results  of  Walsh  and  Donnelly  [8].  Only  one  type  of  boundary  forcing  is  considered;  the  case 
when  the  inner  cylinder  rotates  at  constant  angular  velocity  and  the  outer  cylinder  is  modu- 
lated about  a zero  mean.  In  the  experiments,  Walsh  and  Donnelly  varied  the  gap  width,  the 
forcing  frequency,  and  the  outer  cylinder  angular  velocity  modulation  amplitude.  For  this 
particular  configuration,  it  was  found  in  the  experiments  that  the  modulation  stabilizes  the 
flow  when  compared  to  the  case  with  zero  outer  cylinder  rotation. 

For  the  present  calculations,  only  axisymmetric  disturbances  (n  = 0)  are  considered.  We 
present  the  calculated  results  here  in  terms  of  the  parameters  used  by  Walsh  and  Donnelly 
and  those  used  in  earlier  theoretical  studies  [7].  The  Reynolds  number  Re  used  in  the  previous 
studies  is  related  to  the  Reynolds  number  that  appears  in  the  theoretical  formulation  here 
by  Re  = Re^Jr)/(  1 — 7).  The  dimensionless  forcing  frequency  used  previously,  7,  is  related  to 
the  dimensionless  frequency  of  the  present  formulation  by  u>  = 2y2;  the  use  of  7 rather  than 
u)  compacts  the  frequency  scale.  The  remaining  parameters  7 and  e,  the  radius  ratio  and 
dimensionless  modulation  amplitude,  respectively,  are  defined  in  Section  3.  For  simplicity  in 
all  the  figures  that  follow,  we  have  shifted  the  origin  of  r so  that  r ranges  from  1 to  2 instead 
of  7/(1  - 7)  to  1/(1  - 7). 

In  Figs.  1 and  2,  the  base  flow  as  a function  of  r for  nine  different  fractional  times  covering 
half  a cycle  is  shown  for  two  values  of  modulation  frequency,  7 = 2 and  6,  respectively.  The 
base  flow  velocity  is  normalized  by  the  maximum  velocity  of  the  outer  cylinder.  For  clarity 
only  half  the  cycle  is  shown.  For  the  higher  frequency,  7 = 6,  the  effect  of  the  outer  cylinder 
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modulation  is  confined  to  the  vicinity  of  the  outer  boundary,  whereas  for  the  lower  frequency 
7 = 2,  the  effect  of  the  modulation  spans  the  entire  domain. 

In  the  experiments  of  Walsh  and  Donnelly,  three  radius  ratios  (7  = 0.719,0.88  and  0.95) 
and  two  modulation  amplitudes  (e  = 0.5  and  1.5)  were  investigated  for  several  different  forc- 
ing frequencies,  7.  Calculated  results  are  presented  for  the  same  combinations  of  parameters 
as  the  experiments  for  a larger  and  more  continuous  range  of  frequencies.  In  Fig.  3,  we  show 
the  critical  Reynolds  number  Rc  (the  absolute  minimum  of  Re  as  a function  of  wavenumber. 
a)  versus  7 for  radius  ratio  7 = 0.719  and  modulation  amplitude  e = 1.5.  An  e value  of  1.5 
means  that  the  maximum  amplitude  of  the  periodic  angular  velocity  of  the  outer  cylinder  is 
1.5  times  greater  than  the  steady  angular  velocity  value  of  the  inner  cylinder.  In  the  figure, 
the  curves  represent  the  present  linear  theory  predictions,  while  the  individual  points  are 
from  the  experiments  [8].  Points  below  the  curves  are  stable  according  to  linear  theory,  while 
points  lying  above  are  unstable.  As  can  be  seen,  the  agreement  is  quite  good. 

For  large  frequencies,  the  predicted  critical  values  asymptote  to  the  value  for  a stationary 
outer  cylinder,  since  the  effects  of  the  modulation  are  confined  to  a thin  layer  near  the  outer 
boundary.  The  calculated  curve  approaches  the  critical  value  for  a stationary  outer  cylinder, 
which  for  7 = 0.719  corresponds  to  Rc  = 51.04.  The  theoretical  predictions  consist  of  two 
individual  curves.  The  curve  on  the  right  is  a synchronous  mode,  for  which  the  imaginary 
part  of  the  growth  rate  cr  vanishes.  The  narrow  parabolically  shaped  curve  corresponds  to 
subharmonic  response,  where  <r,-  is  equal  to  w/2.  In  the  experiments,  there  was  no  means 
for  monitoring  the  temporal  character  of  the  instability  as  onset  ensued,  so  the  distinction 
between  synchronous  or  subharmonic  response  could  not  be  made.  The  maximum  amount 
of  stabilization  occurs  for  7 in  the  range  of  3 to  4. 

In  Fig.  4,  we  show  the  critical  values  for  the  case  where  the  radius  ratio  7 = 0.88  and 
e = 0.5.  Again,  the  solid  curve  is  the  linear  theory  predictions  and  the  individual  points 
are  from  Walsh  and  Donnelly’s  experiments.  The  amount  of  stabilization  obtained  in  this 
case  is  significantly  less  than  the  previous  case,  for  which  the  maximum  modulation  angular 
velocity  was  1.5  times  the  constant  angular  velocity  of  the  inner  cylinder.  For  the  range 
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of  frequencies  considered,  only  synchronous  response  is  obtained  at  onset  from  the  theory 
with  this  lower  value  of  e.  The  agreement  between  theory  and  experiment  is  again  quite 
good.  For  7 = 0.88,  the  critical  value  for  the  unmodulated  case  is  Rc  = 44.49,  which  is  the 
asymptotic  value  of  the  theoretical  results  as  the  frequency  becomes  large.  Barenghi  and 
Jones  [9]  present  a curve  of  Rc  versus  7 based  on  their  linear  theory  calculations  for  this  set 
of  parameters.  The  maximum  stabilization  predicted  by  their  calculations  agrees  very  well 
with  our  calculations.  However,  their  calculated  values  of  Rc  do  not  decrease  as  sharply  as 
ours  for  7 greater  than  about  1.25. 

Again  for  7 = 0.88,  Fig.  5 shows  the  onset  conditions  in  the  case  when  the  modulation 
amplitude  is  increased  to  e = 1.5.  The  theoretical  curve  exhibits  similar  structure  as  the  case 
shown  in  Fig.  3,  with  subharmonic  response  again  being  obtained.  Only  five  experimental 
data  points  for  a narrow  range  of  frequencies  were  provided  for  this  case.  At  lower  frequencies 
the  theoretical  stability  boundary  consists  of  parabolic  regions  of  both  subharmonic  and 
synchronous  response,  with  the  overall  trend  being  greater  and  greater  stabilization.  For 
this  case,  the  agreement  between  the  theory  and  experiment  is  less  satisfactory.  Except  for 
one  point,  the  experimental  values  lie  below  the  theoretical  curve.  The  primary  difference 
between  this  case  and  the  one  shown  in  Fig.  3 is  the  frequency  range  of  the  data. 

The  final  comparison  with  the  experimental  data  is  presented  in  Table  1.  For  7 = 0.88 
and  a fixed  frequency  of  7 = 2.3,  the  modulation  amplitude  is  increased  from  e = 0.3  to 
2.0.  For  modulation  amplitudes  below  1.0,  the  agreement  between  theory  and  experiment 
is  good;  however,  for  the  values  e = 1.5  and  2.0,  the  theoretical  predictions  are  considerably 
higher  than  the  experimental  values.  Note  that  the  temporal  response  obtained  from  the 
theory  is  synchronous  for  all  but  the  largest  value  of  e. 

Walsh  and  Donnelly  also  presented  four  data  points  for  7 = 0.95  and  e = 0.5  with  7 in 
the  range  0.6  to  0.9.  The  Rc  from  the  experimental  data  is  approximately  45  within  the 
experimental  error  for  the  range  of  7 covered.  The  theory  predicts  Rc  of  45.8  to  45.9  over 
this  range,  which  is  in  good  agreement  with  the  experimental  data.  The  unmodulated  Rc  is 
42.4  , so  only  a small  amount  of  stabilization  occurs  for  this  case. 
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In  order  to  better  understand  the  results,  and  perhaps  explain  the  discrepancy  for  the 
high  modulation  amplitude  case  shown  in  Fig.  5,  the  structure  of  the  linear  eigenmodes 
was  investigated  both  spatially  and  temporally  for  77  = 0.88,  e = 1.5  and  for  frequencies 
7 = 2 and  6.  For  the  higher  frequency,  Fig.  6 shows  three-dimensional  plots  of  the  u and 
v velocity  components  as  functions  of  the  radial  coordinate,  r,  and  time  (normalized  over 
one  period).  For  this  high-frequency,  synchronous  case,  the  spatial  and  temporal  structure 
of  the  eigenmodes  is  fairly  simple,  having  a single  maximum  in  space  and  being  roughly 
sinusoidal  in  time.  For  7 = 2,  the  critical  Reynolds  number  ( Rc  = 94.0)  corresponds  to  a 
subharmonic  mode,  but  there  is  also  a synchronous  mode  close  by  with  Re  = 95.1.  In  Fig. 
7,  the  velocity  components  for  the  synchronous  mode  are  shown  in  order  to  compare  the 
effect  of  the  frequency  on  the  solution  behavior  for  the  same  type  of  temporal  response.  In 
contrast  to  the  high  frequency  solution,  at  7 = 2,  both  the  spatial  and  temporal  behavior 
are  complicated;  the  temporal  behavior  consists  of  a quiescent  interval  present  for  nearly 
one-half  the  period  followed  by  rapid  increase  in  the  velocity  amplitudes.  In  order  to  better 
view  the  structure,  the  solutions  are  plotted  from  t = 0.25  to  1.25,  so  that  the  nonzero 
regions  appear  in  the  interior  of  the  plots  and  not  at  the  edges. 

Three-dimensional  plots  of  the  velocity  components  for  the  subharmonic  mode  at  7 = 2 
are  shown  in  Fig.  8.  Since  the  temporal  response  is  subharmonic,  the  velocities  are  plotted 
over  two  periods.  Similar  to  the  synchronous  mode  at  the  lower  frequency,  there  are  quiescent 
intervals  followed  by  rapid  change  in  the  amplitudes  over  a short  period  of  time.  In  the 
subharmonic  case,  the  rapid  change  in  amplitude  occurs  once  every  forcing  period,  but 
alternates  in  direction  (positive  or  negative)  so  that  the  solution  is  periodic  over  twice  the 
forcing  period.  For  the  synchronous  case,  the  rapid  change  in  amplitude  is  always  in  the  same 
direction  and  the  period  of  the  flow  is  the  same  as  the  forcing  period.  Here,  the  solutions  are 
plotted  from  t = 0.5  to  2.5  to  better  show  the  behavior.  Given  the  similarity  of  the  temporal 
structure  (except  for  the  change  in  sign  of  the  subharmonic  mode),  it  is  not  surprising  that 
the  Reynolds  numbers  are  approximately  the  same  for  both  modes.  However,  there  is  a shift 
in  the  critical  wavenumber  between  the  two  modes,  with  a = 3.10  for  the  svnchronous  mode 
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and  a = 4.35  for  the  subharmonic  mode. 

Fig.  6 indicates  that  for  large  7,  the  spatial  and  temporal  behavior  is  simple,  whereas  for 
small  7 Figs.  7 and  S indicate  that  the  temporal  behavior  is  complicated.  To  further  illustrate 
this  in  Fig.  9,  we  show  convergence  of  the  temporal  Fourier  coefficients  for  the  v component 
of  velocity.  The  results  correspond  to  the  two  cases  shown  in  Figs.  6 and  7,  representing 
7 = 6 and  2,  respectively.  For  the  larger  7 the  Fourier  coefficients  decay  rapidly,  while  for 
the  smaller  7 a large  number  of  coefficients  are  required  to  obtain  an  accurate  solution.  The 
difficulty  in  obtaining  accurate  solutions  at  low  frequencies  was  also  discussed  by  Barenghi 
and  Jones  [9],  where  they  show  that  Rc  values  below  the  unmodulated  value  can  be  obtained 
with  insufficient  temporal  resolution.  They  attribute  this  as  the  reason  why  the  calculations 
of  Carmi  and  Tustaniwskyj  [7]  predict  destabilization  for  the  specific  boundary  forcing  case 
considered  here. 

While  the  spatial  structure  of  the  eigensolutions  is  less  complicated  than  the  temporal 
structure,  it  is  equally  important  to  assure  that  proper  resolution  is  obtained.  Both  of 
the  numerical  approaches  employed  here  provide  for  high  spatial  accuracy.  The  SUPORT 
approach  achieves  high  accuracy  through  the  use  of  a variable  step-size  high-order  Adams 
integrator.  The  Chebyshev-pseudospectral  representation  of  the  second  solution  approach  is 
also  highly  accurate  as  displayed  in  Fig.  10  by  the  decay  of  the  Chebyshev  coefficients  of 
the  v velocity  for  the  same  7 values  as  Fig.  9.  Typically,  twelve  or  sixteen  Chebyshev  modes 
were  used  for  the  present  calculations,  for  which  the  coefficients  have  decayed  by  six  orders 
of  magnitude.  Clearly,  Fig.  10  also  indicates  that  the  spatial  complexity  does  not  depend 
significantly  on  7. 

The  variation  Re  with  wavenumber,  a,  for  7 = 2,  7 = 0.88,  and  e = 1.5  is  shown  in  Fig. 
11.  As  is  apparent  from  this  figure,  the  minima  of  the  two  lowest  curves  (synchronous  and 
subharmonic)  occur  at  approximately  the  same  Reynolds  number  with  the  absolute  minimum 
on  the  subharmonic  curve.  With  a slight  increase  in  7,  the  absolute  minimum  occurs  on  the 
synchronous  curve  as  is  evident  from  the  transition  in  the  critical  Reynolds  number  plot 
shown  in  Fig.  5.  There  is  also  another  synchronous  branch  at  larger  wavenumbers,  but  the 
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minimum  of  this  branch  occurs  at  a higher  value  of  Re. 

The  effect  of  the  temporal  modulation  on  the  wavenumber  of  the  critical  disturbances 
is  shown  in  Fig.  12,  which  plots  the  critical  wavenumber  as  a function  of  modulation  fre- 
quency, 7,  for  the  same  conditions  as  Fig.  5.  Without  modulation  the  critical  wavenumber 
is  3.12.  For  large  7,  the  critical  wavenumber  asymptotes  to  this  value.  As  7 decreases  the 
critical  wavenumber  decreases  slightly  and  then  rises  sharply.  At  each  of  the  transitions  be- 
tween synchronous  and  subharmonic  response,  there  is  a discontinuous  change  in  the  critical 
wavenumber.  At  small  7,  corresponding  to  large  degrees  of  stabilization,  the  wavenumber 
behavior  is  complicated  with  no  clear  trend. 

7.  Discussion 

Our  numerical  linear  stability  results  are  in  excellent  agreement  with  the  experimental  data 
of  Walsh  and  Donnelly,  except  for  low  frequencies  and  large  modulation  amplitudes  where 
the  experimental  results  he  below  the  calculated  critical  Reynolds  numbers.  Both  the  cal- 
culations and  the  experiments  show  that  modulation  of  the  outer  cylinder  about  zero  mean 
stabilizes  the  flow  when  compared  to  the  case  of  a stationary  outer  cylinder.  For  constant 
rotation  of  the  outer  cylinder,  Taylor  showed  [1]  that  the  critical  Reynolds  number  is  a U- 
shaped  function  of  the  angular  velocity  of  the  outer  cylinder,  and  that  the  lowest  critical 
Reynolds  number  occurs  for  zero  outer  cylinder  rotation.  A simple  argument  for  the  stabi- 
lization by  modulation  of  the  outer  cylinder  assumes  that  the  instantaneous  angular  velocity 
of  the  outer  cylinder  can  be  used  to  characterize  the  instability  [8].  Thus,  modulation  ef- 
fectively moves  the  system  into  a less  unstable  region  of  the  stability  diagram;  in  fact,  for 
sufficiently  large  modulations  the  system  spends  a portion  of  its  time  in  the  stable  region  of 
the  stability  diagram  of  the  unmodulated  flow.  Therefore,  the  stabilization  increases  with 
increasing  modulation  amplitude.  A comparison  of  the  base  flow  (Fig.  1)  with  the  temporal 
response  of  the  eigenfunctions  shown  in  Figs.  7 and  8 indicates  that  the  quiescent  inter- 
val occurs  during  the  half-cycle  when  the  outer  cylinder  velocity  decreases  from  maximum 
velocity  (corotation)  to  minimum  velocity  (counterrotation). 
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There  are  a number  of  possible  explanations  for  the  discrepancy  between  our  calculations 
and  the  experimental  results  for  the  critical  Reynolds  numbers  at  low  frequencies  and  large 
modulation  amplitudes  (see  Fig.  5 and  Table  1).  Walsh  and  Donnelly  indicate  that  at  low 
frequency  it  becomes  increasingly  difficult  to  determine  the  onset  of  instability.  In  linear 
stability  calculations,  the  criterion  for  stability  is  that  of  transient  stability  [IS]  where  a 
disturbance  may  grow  for  part  of  the  forcing  cycle  but  ultimately  decays.  At  the  onset  of 
transient  instability,  the  eigenmodes  (Figs.  7 and  8)  show  a quiescent  interval  followed  by 
a very  large  amplification.  Depending  on  the  size  of  initial  fluctuations,  this  large  amplifi- 
cation might  invalidate  the  linear  theory,  or  alternatively  might  be  viewed  as  instability  in 
experimental  observations  ([18],  [19]).  Since  disturbances  are  always  present  in  experiments 
and  since  individual  disturbances  may  grow  to  significant  amplitude  for  part  of  the  cycle, 
this  situation  may  be  viewed  as  a manifestation  of  secondary  flow,  even  though  an  individual 
disturbance  would  decay  over  several  cycles  according  to  linear  theory.  The  likelihood  of  this 
situation  occurring  is  greater  at  lower  frequencies,  where  the  time  period  for  the  effect  of 
multiple  disturbances  to  accumulate  is  longer. 

Finally,  it  is  important  to  mention  that  although  the  present  study  is  comprehensive,  it 
is  not  exhaustive  owing  to  the  large  number  of  parameters  and  the  complexity  exhibited  by 
the  stability  boundaries.  For  example,  it  was  not  possible  to  evaluate  Re  for  all  wavenumbers 
in  the  entire  frequency  range.  As  shown  in  Fig.  11,  there  are  multiple  minimums  in  the  plot 
of  Re  versus  a at  fixed  frequency.  Also,  only  axisymmetric  disturbances  were  investigated 
here.  It  is  possible  that  for  certain  parameter  values  nonaxisymmetric  modes  may  be  more 
unstable.  It  is  known  that  for  steady  rotation  the  critical  modes  can  be  nonaxisymmetric 
when  the  cylinders  rotate  in  opposite  directions  [20],  However,  given  the  good  agreement 
with  the  experimental  results  for  all  but  low  frequency,  high  amplitude  modulation,  it  is 
clear  that  axisymmetric  linear  theory  is  still  relevant  for  studying  the  effect  of  time-periodic 
forcing  in  the  Taylor-Couette  problem. 
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Table  I.  Comparison  of  Rc  from  the  experimental  data  of  Walsh  and  Donnelly  and  the  present 
calculations  for  various  modulation  amplitudes,  e,  at  fixed  frequency,  7 = 2.3,  with  7 = 0.88. 
The  calculated  wavenumber,  a,  and  mode  type  are  also  shown. 
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0.0 

0.3 
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Rc  (Data) 

44.5 
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Rc  (Calc) 

44.5 
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Figure  Captions 

Figure  1.  Azimuthal  base  flow  velocity,  u(0)  as  a function  of  the  radial  coordinate,  r,  for 
various  fractional  times  through  half  a period  for  7 = 2,  7 = 0.88,  and  e = 1.5. 

Figure  2.  Azimuthal  base  flow  velocity,  v^0)  as  a function  of  the  radial  coordinate,  r,  for 
various  fractional  times  through  half  a period  for  7 = 6,  7 = 0.88,  and  e = 1.5. 

Figure  3.  Critical  values  of  Re  as  a function  of  7 for  7 = 0.719  and  e = 1.5.  The  points 
are  the  experimental  results  of  Walsh  and  Donnelly;  the  curves  are  the  present  linear  theory 
(solid  curve  - synchronous  response,  dotted  curve  - subharmonic  response). 

Figure  4.  Critical  values  of  Re  as  a function  of  7 for  7 = 0.88  and  e = 0.5.  The  points  are 
the  experimental  results  of  Walsh  and  Donnelly;  the  solid  curve  is  the  present  linear  theory. 

Figure  5.  Critical  values  of  as  a function  of  7 for  7 = 0.88  and  e = 1.5.  The  points 
are  the  experimental  results  of  Walsh  and  Donnelly;  the  curves  are  the  present  linear  theory 
(solid  curve  - synchronous  response,  dotted  curve  - subharmonic  response). 

Figure  6.  Linear  eigenfunction  velocity  components  u(r,  t)  and  v(r,  t)  plotted  for  one  forcing 
period  over  the  spatial  domain  for  7 = 6,  7 = 0.88,  and  e = 1.5  (synchronous  response). 

Figure  7.  Linear  eigenfunction  velocity  components  u(r , t)  and  u(r,  t)  plotted  for  one  forcing 
period  (starting  at  t = 0.25)  over  the  spatial  domain  for  7 = 2,  7 = 0.88,  and  e = 1.5 
(synchronous  response). 

Figure  8.  Linear  eigenfunction  velocity  components  u(r,  t)  and  u(r,  t)  plotted  for  two  forcing 
periods  (starting  at  t = 0.50)  over  the  spatial  domain  for  7 = 2,  7 = 0.88,  and  e = 1.5 
(subharmonic  response). 

Figure  9.  The  m th  temporal  Fourier  coefficient  of  the  v velocity  component  at  a fixed 
spatial  point  as  a function  of  m for  the  linear  eigenfunctions  shown  in  Fig.  6 (solid  curve) 
and  Fig.  7 (dotted  curve). 
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Figure  10.  The  nth  Chebyshev  spectral  coefficient  of  the  v velocity  component  at  a fixed 
time  as  a function  of  n for  the  linear  eigenfunctions  shown  in  Fig.  6 (solid  curve)  and  Fig. 
8 (dotted  curve). 

Figure  11.  Values  of  Re  as  a function  of  wavenumber,  a,  for  7 = 2.0,  7 = 0.88  and  e = 1.5. 
(solid  curve  - synchronous  response,  dotted  curve  - subharmonic  response). 

Figure  12.  Critical  values  of  a as  a function  of  7 for  the  same  parameters  as  Fig.  5.  (solid 
curve  - synchronous  response,  dotted  curve  - subharmonic  response). 
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